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ABSTRACT 

A numerical method is developed for solving periodic, three-dimensional, vortical 
flows around lifting airfoils in subsonic flow. The first-order method that is presented fully 
accounts for the distortion effects of the nonuniform mean flow on the convected upstream 
vortical disturbances. The unsteady velocity is split into a vortical component which is 
a known function of the upstream flow conditions and the Lagrangian coordinates of the 
mean flow, and an irrotational field whose potential satisfies a nonconstant-coefficient, 
inhomogeneous, convective wave equation. Using an elliptic coordinate transformation, 
the unsteady boundary value problem is solved in the frequency domain on grids which 
are determined as a function of the Mach number and reduced frequency. Extensive 
comparisons are made with known solutions to unsteady vortical flow problems, and it is 
seen that the agreement is in general very good for reduced frequencies ranging from 0 up 
to 4. 



I. INTRODUCTION 


Most flows encountered in aerodynamics are high speed flows where the Reynolds 
number is large and the effects of viscosity are confined to small regions such as boundary 
layers and wakes. Because major portions of these flow fields are essentially inviscid and 
irrotational, potential flow theory has been used extensively by aerodynamicists in the 
analysis of flows about streamlined bodies. Today steady potential flow solvers are widely 
used in the design of aircraft wings, turbomachinery blades, and helicopter rotors. 

In many real flow applications, however, the flow is not steady but unsteady. Fre- 
quently the unsteadiness in the flow is due to the occurrence of upstream vortical distur- 
bances that are convected downstream and induce an unsteady flow field as they interact 
with the body. For an aircraft wing, such upstream flow distortion can be caused by at- 
mospheric turbulence. For propeller and turbomachinery blades, the vortical disturbances 
may be caused by the viscous wakes of an upstream rotor or stator, installation effects, or 
upstream turbulence. 

When viewed from the blade frame of reference, the upstream vortical disturbances 
will appear as propagating vorticity waves that are called gusts. There are a number of 
undesirable effects that can be associated with such vortical gusts. They will, for example, 
induce unsteady forces on the airfoil surface which can cause forced vibrations and radiate 
noise into the far field. In some instances, the impinging gusts may cause flow separation 
and loss of aerodynamic performance. For rotating blades, the fundamental frequency of 
the upstream disturbances will equal the blade passing frequency. If the frequency of the 
aerodynamic excitation equals a natural frequency of the rotating blades and the amplitude 
is sufficient, then catastrophic structural failure may result. 

Another possible source of unsteadiness in the flow is the unsteady motion of the 
airfoils or blades themselves. Such unsteady structural motion can be caused by structure- 
borne vibrations as well as the flow-induced oscillations described above. There can also 
be unsteady interactions between the airfoil motion and the incident disturbances which 
can dampen or increase the magnitude of the airfoil unsteady motion. 

Because of the undesirable effects associated with these unsteady flows, there is con- 
siderable interest in controlling and understanding the aerodynamic excitations which can 
cause such unsteady blade motion. 

Up until recently, most numerical efforts to solve these kinds of unsteady flows con- 
centrated on potential methods. The early work dealt with solving the unsteady small 
disturbance potential equation as a way of obtaining the unsteady flow around oscillating 
airfoils or cascades. Later work was directed toward solving the linearized unsteady po- 
tential equation and the unsteady full potential equation. Potential methods have proven 
to work well for oscillating airfoil problems, but unfortunately they cannot adequately 
account for the vortical part of the flow. Previous potential formulations which have 
included the effects of the upstream vorticity have invoked the linear thin airfoil approxi- 
mation and assumed that the imposed vortical gust is convected without distortion by the 
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nonuniform mean flow. This was the approach used by McCroskey and Goorjian 1 and 
McCroskey 2 . However, as shown by Goldstein and Atassi 3 , Atassi 4 , and Scott and 
Atassi 5 , the assumption that the gust is convected without distortion is not justified and 
is a poor approximation for flows with a spatially varying mean flow. This is especially 
true for turbomachinery and propeller flow fields where the blades are heavily loaded and 
there are strong mean flow gradients. 

In the past few years, computational efforts in unsteady aerodynamics have concen- 
trated on the so-called primitive variable methods in which the unsteady Euler or Navier- 
Stokes equations are solved in time along with certain specified boundary conditions. Un- 
like the potential methods, the primitive variable methods are equally well-suited to both 
oscillating airfoil problems and flows with convected upstream vorticity. The main diffi- 
culty associated with these methods is that they are too expensive to be used for routine 
engineering calculations such as are encountered in design work. In addition, uncertain- 
ties about physically correct far field boundary conditions leaves some question as to the 
accuracy of the solutions. 

In a previous paper 5 , the authors presented a linearized unsteady aerodynamic anal- 
ysis for subsonic vortical flows around lifting airfoils which represents an alternative to the 
potential and primitive variable methods. The analysis that was presented in [5] offers 
the computational efficiency of potential methods, but at the same time accounts for the 
convection and distortion of the upstream vorticity by the nonuniform mean flow. Our 
analysis is therefore equally well-suited to vortical flow problems as well as to oscillating 
airfoil problems. In addition, since our linearization is about the nonuniform mean flow, 
the full nonlinear effects of the mean flow are accounted for. Only the unsteady part of 
the flow is linearized. These features, coupled with the inherent efficiency of the linearized 
approach, make our analysis an ideal solution method for unsteady aerodynamic flow fields. 

In [5] we presented the mathematical formulation of the general linearized boundary 
value problem, and demonstrated the capabilities of our approach by presenting numerical 
solutions for a large variety of unsteady vortical flows. The numerical results presented 
showed in detail the effects of airfoil thickness, angle of attack, camber, and Mach number 
on the unsteady lift and moment of isolated airfoils subjected to periodic vortical gusts in 
subsonic flows. 

In the present paper our major purpose is to present the details of the frequency 
domain numerical scheme that has been developed to implement our linearized unsteady 
aerodynamic analysis. Since we presented the full formulation of the analysis in [5], we 
will only summarize the resulting boundary value problem in the present paper. This will 
be done in the following section. Following that will be the main section of the paper 
which presents our frequency domain numerical scheme. Finally, in Section IV, we present 
a large variety of numerical results which demonstrate the validation of our codes. 
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II. LINEARIZED UNSTEADY AERODYNAMIC FORMULATION 


Consider an inviscid, compressible flow past a two-dimensional airfoil placed at 
nonzero incidence to a stream with uniform upstream velocity Uoo in the x\ direction. 
We assume in the present discussion that the fluid is an ideal, non-heat conducting gas 
with constant specific heats, and that there are no shocks in the flow. Under the above 
assumptions there will be a steady potential flow around the airfoil so that we may write 

U 0 (£) = V$o, (2.1) 


where 0 subscripts are used to denote steady mean flow quantities. 

Let us assume that far upstream an unsteady vortical disturbance is imposed on the 
flow. The only restriction that we place on the upstream disturbance is that it can be 
expressed as a generalized Fourier integral so that we may write 


u e 


,(X Z U r~ 


>*) = L 

Jk 


a(k)e ihl - l: - M ~'Zdk 


(2.2) 


and in addition, that its length scale V and characteristic velocity are such that the 
condition 

< — (2.3) 


U 0 


r_ 

Uoo 


is satisfied, where c is the airfoil chord length. We thus require that the time scale associ- 
ated with the mean flow be an order of magnitude less than the time scale associated with 
the upstream unsteady disturbances. Since our concern is with flows that have large scale 
upstream disturbances in which l' is the same order of magnitude as the chord length c, 
condition (2.3) essentially reduces to the requirement that Mqo <C Uoo- 

Since we present a linearized mathematical formulation, we may without loss of gen- 
erality consider a single Fourier component of the incident vortical disturbance, and solve 
for more general disturbances by superposition. We therefore consider incident vortical 
gusts of the form 


u c 


_ fciWX-tUcot) 


(2.4) 


where a and k must satisfy 


a ■ k = 0 


(2.5) 


to ensure that the continuity equation is satisfied. 

In general, the components of X in equation (2.4) are not the spatial coordinates, but 
rather are Lagrangian coordinates of the mean flow. For the case of two-dimensional mean 
flow the components of X are given by 


X 2 = 


PoqU oo 


(2.6) 
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and 


•^3 = (2-7) 

where is the stream function of the mean flow and x 3 is the spatial coordinate in the 
spanwise direction. The component X\ is defined by 


X x = UooA, 


( 2 . 8 ) 


where A is the Lighthill 
as 


“drift” function 6 , which can be expressed in terms of $ 0 and T 0 



where the integration is carried out on tHo = constant. The difference in A between two 
points on a streamline is the time it takes a mean flow fluid particle to traverse the distance 
between those two points. Note that for the thin airfoil case in which the mean flow is a 
uniform parallel flow, the components of X reduce precisely to the spatial coordinates. 

We assume that the total unsteady flow field can be represented by 


U(x,t) = Uo(x) + u(x,t ) 

(2.10) 

p(x,t) =po(x)+p'(x,t) 

(2.11) 

p(x, t ) = p 0 (x) + p'(x, t) 

(2.12) 

s(x, t) — Sq + s'(x, t) 

(2.13) 


where the entropy s 0 is constant, and u, p' , p' , and s' are the unsteady perturbation 
velocity, pressure, density and entropy, respectively. Quantities with 0 subscripts are the 
steady mean flow quantities which are assumed to be known. Note that these quantities 
obey the steady nonlinear equations of motion, so that the linearization of the unsteady 
flow is about the fully nonlinear mean flow. 

Substituting relations (2.10) - (2.13) into the nonlinear Euler equations and neglecting 
products of small quantities, one obtains the linearized continuity, momentum, and energy 
equations 

+ p'V ■ U 0 + V • ( Po u) = 0 (2.14) 

Po(^ + u- Wo) + p'Uo • Wo = Vp' (2.15) 

¥=°’ (2 - i6) 
where ^ = A 4- Uo ■ V is the convective derivative associated with the mean flow. 

Dt at u 

It is shown in References 7 and 8 that if the unsteady velocity is decomposed into 
the sum of a known vortical component and an unknown potential component, then the 
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problem for determining the unsteady flow may be reduced to solving a single, non-constant 
coefficient, inhomogeneous convective wave equation which may be written 


Do. 1 Do(f>. 


-V.(p 0 V < A)=-V.(p 0 ^ ) ), 

Po po 


(2.17) 


where 

u(x, t) = u^ RS> + V^>. 

The imsteady pressure is related to (f> through the equation 


p = -po(x) 


Do<t> 
Dt ' 


(2.18) 


(2.19) 


The vortical velocity u <R> is a known function of the mean flow Lagrangian coordinates 
and the upstream vortical disturbances and is given by 


u 


(R) 


[V(3-A')]e 


ik-(X-iUooi) 


+ V«^, 


( 2 . 20 ) 


where 


<f> 


■(oi + 


q 2 fci - aik 2 1 - e lkA yk-jx-iu^t) 
1 + iaoUooki h 2 


( 2 . 21 ) 


and 


a = («i,a2) a 3) an d 


,dU 0 - 1 

do = -(-w— ) • 

on s 


( 2 . 22 ) 


Here n denotes the direction of the outward unit normal, S denotes the stagnation point 
near the airfoil leading edge, and Uq = |^o| is the magnitude of the mean velocity. (See 
[5] or [8] for more details concerning the purpose and derivation of the function 4> ■) 
Finally, the potential <f) must satisfy the boundary conditions 


V(f> • n = 0 

airfoil surface 

(2.23) 

IS< A « = ° 

wake 

(2.24a) 

A[V<^> • n] = 0 

wake 

(2.24b) 

V<f) — > — V<^> as 

Xi — > — oo, 

(2.25) 


where equation (2.23) is the impermeability condition at the airfoil surface, equations 
(2.24a) and (2.24b) impose continuity of the pressure and normal velocity across the wake, 
and equation (2.25) ensures that u(x,t ) — > u^x^t) as xi — > — oo. 

The linearized boundary value problem for the unsteady gust response problem thus 
consists of the governing equation (2.17) and boundary conditions (2.23) — (2.25), together 
with the requirement that (f) is continous at the airfoil trailing edge. 
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III. NUMERICAL SCHEME 


Reformulation and Nondimensionalization of the Boundary Value Problem 

For numerical purposes it is necessary to reformulate the boundary value problem 
presented in the previous section into a form more suitable for numerical computations. 
Of particular concern is condition (2.25). In order to facilitate the implementation of the 
far field boundary condition, it is convenient to replace <j> by a function whose gradient 
vanishes as r — » oo, where r is the distance from the airfoil center. 

To this end, we introduce the potential functions <f>i and <j> 2 , where 

$ = <}> i-fa (3.1) 

and 4>2 is a known function which is constructed such that 

\<j >2 — — > 0 as r — > 00 . (3-2) 

Equation (3.1) together with conditions (2.25) and (3.2) then show that the new potential 
function (j) 1 will satisfy 


V^i — » V (j >2 — V</> —■ y 0 as r — > 00 . 


(3.3) 


The problem may then be reformulated in terms of the unknown potential <f>\. 
To satisfy condition (3.2), the function (j> 2 must take the form 


^ + 1 -HaoCJt, 


a 2 k x -a x k 2 1 — e lk * x * , 


-)< 


(3.4) 


where the vector X e satisfies 


\X e — X\ — > 0 as r — > 00 . 


(3.5) 


To satisfy this condition for the general problem of vortical flows around lifting airfoils, we 
define X e as follows: 


X, 


e,l 


$0 

ul 


r 7T 

^-sgn(T 0 )[- +sgn('k 0 )tan' 






/loo'k 


0 


] (3.6) 


X e ,2 = X 2 (3.7) 

X e ,3=X 3 . (3.8) 

The expression for X e> i is obtained by making a far field expansion of Lighthill’s drift 
function A in terms of $0 an d 4V Note that the first term in the expansion is just 
and that the second term arises due to the circulation around the airfoil. Since the second 
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term vanishes for airfoils with zero circulation, it is clear that the formulation of the source 
term for nonlifting airfoils is much simpler than for lifting airfoils. 

We also point out that the first factor in brackets in equation (3.6) is discontinuous 
and undefined at the points on the airfoil where T 0 = 0,T 0 = 0 + and $ 0 = 0,T 0 = 0~. 
The second factor in brackets is not part of the expansion itself, but is included to remove 
the discontinuity in X e p. By including the second factor in brackets and defining X e ,i = 0 
at $ 0 = 0, To = 0, we obtain an expression for X e p which is everywhere continuous. 
It is important that X e p be continuous along the airfoil surface, for if it were not, the 
potential function <j>\ would have to satisfy a discontinuous airfoil boundary condition. 
(See condition (3.10) below.) By defining X e p as in (3.6), we obtain an expression for X e 
which is everywhere continuous and also ensures that conditions (3.5) and (3.2) will be 
satisfied for both lifting and nonlift ing airfoils . 

Finally, the parameter /?oo = \/l — Moo 2 , where M 0 0 is the free stream Mach number, 
and c is the airfoil chord length. The term arises due to a Gothert’s rule correction 
on the mean velocity so that the expression for X e p is valid for both compressible and 
incompressible flows. 

Before presenting the reformulated boundary value problem in terms of the potential 
(f>i , we present the nondimensionalization of the problem. We normalize as follows: 


*1,^2, *3, Xi,X2,Xs, X e ,l 

by 

c 

2 

$ 0 ,r 

by 

-U 

2 u oo 

To 

by 

2 pooU oo 

Uo, co 

by 

tfoo 

Po 

by 

Poo 

P' 

by 

Poo^oo |®| 

t, A 

by 

c 

■zu^ 

u> 

by 

2 U oo 
c 

k\ i k 2 , &3 

by 

2 

c 


by 

m 

— * 
a 

by 

|a| 


The normalized wave number k\ = ^rr - , where u and are the dimensional angular 
frequency and free stream velocity, respectively, is called the reduced frequency. 

We will assume throughout the remainder of the present section, unless stated other- 
wise, that all quantities are nondimensional. 

The governing equation for <j > i is then 

> - "V ■ (p„V« = W ■ (p„u™) 

Dt c 0 z Dt p o p o 

+ ^(-%^)--W-UW 2 ) (3.9) 

Dt Co Dt p o 
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and the boundary conditions are 


Do 

Dt 


V<^i • n = 'S/<f >2 • n 

[A(^! - <f> 2 )] = 0 


airfoil surface 


wake 


A[V(</>i — <j> 2 ) • n\ = 0 wake 

W(j) 1 — > 0 as xi — > — 00 . 


(3.10) 

(3.11a) 

(3.11b) 

(3.12) 


For completeness we also present the nondimensional expressions for the potential 
functions and (f> 2 , for the unsteady velocity and pressure, and for the upstream velocity 
disturbances. 


7 */ . a 2 k x -a x k 2 l-e 

* k x ^ 1 + 1 + iaoh k 2 


ik-X —ik\ t 


(3.13) 


where 


a 2 ki -aik 2 l-e yk-x.-ik^t^ 


h -h {ai+ 1 +«.„*, k 2 


A', — A, X 2 — X «, 2 — $o. A 3 = A ,, 3 = X 3 

r ,.7r . $0 


(3.14) 


(3.15) 


Xe,\ = $0 - — sgn(^o)["j + sgn(T 0 )tan ^ )][1 - e (<i>0+ ' 1 '° ) ] (3.16) 

n ^ Poo * 0 


where 


u(x,t) = rf R '> +V(<f> 1 -<I> 2 ) 
u (R) = [V(a • X)]e* tjl - iklt + V4>. 


(3.17) 

(3.18) 


, /-.s D 0 ((j>i — (j> 2 ) 

P = -P.(x) m . 


_ £ ik-X — ik\t 


Urv, = ae 


(3.19) 

(3.20) 


Determination of Mean Potential Flow 

In order to obtain numerical solutions to equation (3.9) and its associated boundary 
conditions, one must first obtain the steady potential flow about the airfoil for the given 
flow conditions. This will in general require the use of a standard potential flow solver 
such as FL036. 9 

However, an examination of equations (3.13) through (3.18) indicates that the most 
natural choice of independent variables in which to solve equation (3.9) are T 0 and To, 
the mean flow potential and stream functions. Since standard potential flow codes solve 
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the steady problem in terms of the spatial coordinates x\ and X 2 , there is some difficulty 
in obtaining the steady solution as a function of To and To. 

Another difficulty arises due to the fact that the grids used by steady flow solvers are 
not suitable for the unsteady calculation. As reported in References 10 and 11 , accurate 
solution of equation (3.9) over a large range of flow conditions requires using grids which 
are determined as a function of both the reduced frequency k\ and the free stream Mach 
number M 0 Q . This means that in general it will be necessary to interpolate the solution 
from the steady grid onto the appropriate unsteady grid. 

Because of the loss of accuracy that can result from such an interpolation process, 
and also because of the need to know the mean flow as a function of T 0 and To, an 
analytical scheme that can obtain the compressible, subsonic flow about isolated airfoils 
was developed. The scheme is based on the idea that, except for a small inner region 
surrounding the airfoil, the flow gradients are not too large. Thus in the large outer region 
extending to infinity, the mean flow is essentially governed by a set of linear equations. As 
a result, one can use Gothert’s Rule, whereby the compressible flow about a given airfoil 
can be obtained from the incompressible flow about a similar airfoil. 

If we let aci @c, and 7 c denote the angle of attack, thickness ratio, and camber ratio of 
the given airfoil in a compressible flow, then the transformed airfoil for the incompressible 
flow field has angle of attack, thickness ratio, and camber ratio given by 


07 = Poo 01 c 1 

0i = Poo0c > (3.21) 

7 / = Poo 7C ) 

where I subscripts denote quantities from the incompressible flow field. Using dimensional 
quantities and denoting the compressible velocity by (U^o +«c, v c) at the point (x, y ) and 
the incompressible velocity by ({Too + u/,17) at the point (37, y/), the spatial coordinates 
and velocity in the compressible and incompressible planes are related by 


X = Xj 

(3.22a) 

yi 

y = J~ 

Hoo 

(3.22b) 

ui 

Uc= H * 

r OO 

(3.23a) 

VI 

vc = ~r- 

Hco 

(3.23b) 


It is assumed here that the free stream velocity U 0 0 is aligned with the x axis, and that 
the angle of attack, thickness ratio, and camber ratio of the airfoil are such that the 
perturbation velocities are small compared to U 0 c . The potential and stream functions of 
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the compressible flow field are then related to the potential and stream functions of the 
incompressible flow field by 


$o - UoqX = — ($/ - UooXi ) (3.24a) 

f'OO 

^0 - - ^ooy/) (3.24b) 

Poo 

Using (3.22), and assuming that all quantities have been nondimensionalized as in the 
previous section, equations (3.24) can be rewritten 

- M^xr = /?^$ o (3.25a) 

^i = /? oo^o (3.25b) 

The problem is then, given (^oj^o)? solve equations (3.25) for ($/, 4//) and then 
use relations (3.22) and (3.23) to determine the spatial coordinates (x,y) and velocity 
components (uc,vc) of the compressible flow field. If this can be done, then we have the 
compressible flow field determined as a function of ($ 0 , 4 r o)- Note that this assumes that 
we can determine (xj,yj) and (uj,Vf) as functions of ($j,Tj). Since there is a one to 
one correspondence between (xj, iji) and ($/, Tj), and between (uj, vj) and ($/, 4/j), this 
is not a difficulty. At the least, it can always be done numerically. For the special case 
of Joukowski airfoils, however, it is possible to express the complex potential (<&j, T/) in 
terms of the polar coordinates (r, 6), and then obtain (a?/,y/) and ( ui,vj ) through known 
functional expressions of (r, 0). 

For the case of incompressible flow around a Joukowski airfoil the complex potential 
is given by 


2 ia/ p /■! 

+ i*i = Cz~' ai + ln( a ) + K (3 ‘ 26) 

where K is an arbitrary constant and 

C' = C - Co' = re ie . (3.27) 

Here (o' = — e + ie' is a complex constant, and the parameters a, e, and e' depend on the 
airfoil geometry. Tj is the steady circulation around the airfoil and is given by 


where /? is defined by 


Tj = — 47rasin(o;/ + /?) 


sin/? = 


a 


(3.28) 


( 3 . 29 ) 
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The spatial coordinates (x /, y / ) are determined from r and 0 through the Joukowski trans- 
formation 

d 2 

xi + i yi = C + ~r (3.30) 


where the parameter d satisfies 


(e + d) 2 + e' 2 = a 


(3.31) 


Finally, the velocity components («/,«/) are given by 

C 2 [(' + ae ! '0' + ^] 


uj — i vi 


(C + d)r 


- 1 . 


(3.32) 


Using relations (3.26) through (3.30), equations (3.25) can be expressed in terms of r 
and 0 as 


(r + ~r) cos(0 — ai) — [2a sin(aj + j3))0 + I\ 

-MUr cos 0- e + {rcos g d l[yl\t7ne+e'y ] = Plo*o t 

(r - ^)sin(0 - ai) + 2asin(o/ + ^)ln(^) = /Joo’f'o , 

If for given $ 0 and T 0 we can solve equations (3.33) for r and 0 , then equations (3.30) 
and (3.32) can be used to get ( xi,yi ) and (u/,v/), and equations (3.22) and (3.23) can 
be used to obtain ( x,y ) and ( u c ,v c ). Once we have obtained (u c ,v c ), the other mean flow 
quantities can then be obtained from Bernoulli’s law for polytropic gases. 

We note that, while the system of equations (3.33) is highly nonlinear in the unknowns 
r and 0, it can be routinely solved by two-dimensional Newton iteration. Once a subroutine 
has been developed to solve equations (3.33), the compressible steady flow around any 
Joukowski airfoil can be very efficiently obtained. In addition, the mean flow is obtained 
for arbitrary $o and To, so that there is no restriction whatsoever on the particular grid 
that may be used for the unsteady calculation. 

The only limitation in obtaining the mean potential flow by this particular approach 
is the underlying assumption that uc and vc must be small compared to Uoo. This means 
that the method will not give a good approximation in the inner region and particularly 
near the stagnation point where the perturbation velocities are of the same order of mag- 
nitude as Uoo- However, extensive testing of this particular approach and comparing with 
the steady potential flow solver FL036 has shown that the region of inaccuracy is very 
small. Figures 1 through 4 show Mach number comparisons between the present approx- 
imate analytical scheme and FL036. The comparison is made at grid points along fixed 
grid lines used by FL036. It is seen that the agreement overall is quite good, with the 



66 



exception of grid points on the airfoil surface that are near the stagnation point. Because 
of this inaccuracy, we use FL036 to calculate the mean flow quantities along the airfoil 
surface itself, and use the approximate analytical scheme off the airfoil except in a small 
region just upstream of the stagnation point. In this region, for airfoils that have steady 
loading, the velocities are calculated using a Taylor series expansion. For airfoils with- 
out steady loading, the velocities are calculated from a local analytical solution which is 
patched to the outer solution. 

Frequency Domain Formulation 

An inspection of equations (3.13), (3.14), (3.18), and (3.20) indicates that the time 
dependence of the present boundary value problem comes entirely through the harmonic 
term e~ tklt . It is therefore possible to make a transformation from the time domain into 
the frequency domain by a simple change of dependent variable. By transforming the 
problem into the frequency domain, time is completely eliminated from the problem and 
it is possible to significantly simplify the mathematical formulation of the boundary value 
problem. 

For the case of two-dimensional mean flow, we transform into the frequency domain 
by making the following change of dependent variable: 

<f> x = <pe~ iklt+ik 313 (3.34) 


By including the i k :i x 3 term in the transformation, the harmonic dependence on the span- 
wise component x% is also eliminated, since all of the e lk3Xa terms then factor out from 
each side of the equation. This is of course possible in view of (3.20) and (3.15). 

Before presenting the governing equation in the new dependent variable ip, we intro- 
duce the linear operators £ and £ 0 to simplify the notation. 


and 


where 


D 0 1 D 0 1 -> -> 

Dt c 0 2 Dt p 0 V ’ ^ 

(3.35) 

„ ,,2 D l , 9 2 d 2 d 2 , 

Co ~ 00 Dt 2 + dtfg + dx\ ^ 

(3.36) 

Do d d 

Dt 0 ~ dt + d$o 

(3.37) 


The governing equation then takes the form 


£<I>1 = — V ■ (p 0 rt R) ) + £fa 

Po 


(3.38) 


The operator £q is essentially the operator for the linear thin airfoil gust response problem. 
By writing the governing equation in the following equivalent form, the left hand side of 
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the equation will exhibit the character of the thin airfoil operator in the far field since 
C — > £o there, 

£a<f>i + (£- £ 0 )<h = —V • ( Po ^ R) ) + Cfo (3.39) 

P o 

We then have in terms of ip 


£o<Pi + (•£ ~ Co)<h — 


+ |V + 2iki Mi ^ - * 2 v] 


r f2 w2 k l M 2 ■, r r2 ^ ,M 2 A1 

+ [ fc l M oo f/2 * fc l ^0 a$Q ( V 2 )lV> 


(3.40) 
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a 2 (p £/ 0 2 dpo dip 
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where M is the local Mach number of the mean flow. 
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Equation (3.40) may be simplified further by making the following change of both 


dependent and independent variables: 


<p = ip e ~ iKo * 

(3.41a) 

where 

A °= p 
r oo 

(3.41b) 

and 


$ = $ 0 

(3.42a) 

^ = ^oo^O 

(3.42b) 


Expressing equation (3.40) in terms of ip and the new independent variables $ and 'I', 
one gets 

£oPi 4 - (£ - £q)4>\ = 


e -ikit+ik s xs e ~i 
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(3.43) 
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where A 1 ...A 5 are functions of ($, T) defined by 


Ai($, \E r ) = 


PL 


k\M 2 

~W 


tt2 d 

lklU °d^u^ 


(3.44a) 
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(3.44b) 


A 3 (*,tt) 


Uldpp_ 

Po <9<I> 


(3.44c) 


a 4 ($, *) = ffl - ful 


(3.44d) 




(3.44e) 


Note that in the far field the functions A 1 ...A 5 tend to zero so that the right hand 
side of (3.43) reduces to that of the linear thin airfoil theory. 12 

To complete the frequency domain formulation, it is necessary to present the right 
hand side of the governing equation (3.39). We will refer to the right hand side of (3.39) 
as the source term and denote it by e.~ tklt+zk 3 Xa S. To facilitate the presentation of the 
source term, we will write the governing equation as 


£0 <h+(C-Co)<h= e" ifclt+ifc »*» S 


(3.45a) 


or 


Co <h+(£- Co) <h = e- iklt+ik * x * (S 1 + S 2 + S 3 - S 4 ) (3.45b) 
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where 


g—ikit+ik 3 X3 c> 


V Po 
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(3.46a) 
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(3.46b) 
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The frequency domain governing equation is then given by 
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(3.46d) 


a?/> , a?/> . a 2 v> , a 2 t/» 

+ ^ + ^2— + a 3 — + a 4 ^ + a 5 ^ 


(3.47) 


= e^o^Si + S 2 + S 3 - S 4 ) 


where Ai ...A 5 are defined by equations (3.44) and Si-..S 4 are defined by equations (3.46). 
In the far field both the coefficients Ai...A 5 and the source term Si + S 2 + S 3 — S 4 tend 
to zero so that the equation reduces to a Helmholtz equation. 

We conclude this section by presenting the frequency domain formulation of the airfoil 
and wake boundary conditions given by equations (3.10), (3.11a.) and (3.11b). In terms 
of the new coordinates ($, T) and the new dependent variable ?/>, equation (3.10) for the 
airfoil boundary condition becomes 


dif) 


JK 0$ 


Qi^ooT 1 — e ® a 2 + ia^aik^ , ik t x el 


Poo 7T $ 

The wake boundary condition (3.11a) becomes 


1 + iaoki ' 


(3.48) 


(-H 1 + U 0 2 ^)lA»e- iK °* - vj )]= 0 , 


(3.49a) 


where 
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(3.49b) 
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and where A(?/>e lK °® — </? 2 ) denotes the jump in the quantity (ipe tK °® — ¥ 2 ) across 
the vortex sheet behind the airfoil. Finally, condition (3.11b) becomes 

AfV^e-^ 0 * - <p 2 )-n\ = 0. (3.49c) 


Transformation into Computational Coordinates and 
Formulation of the Numerical Boundary Value Problem 

Our basic numerical approach to solving equation (3.47) is to use the method of finite 
difference approximations. By discretizing the flow field and employing finite differences 
at each grid point, a large linear system of equations is obtained which can be solved using 
a matrix solver. 

Previous experience in solving equation (3.47) for the case of flat plate and symmetric 
airfoils has shown that the independent variables (4>, T) are not suitable computational 
coordinates for the gust response problem . 10,11 There are difficulties in obtaining consis- 
tently accurate results over a large range of Mach numbers and reduced frequencies, and 
also problems with the implementation of far field boundary conditions. A transformation 
of the independent variables is needed which not only provides an adequate distribution of 
grid points around the airfoil in the near field, but also provides a distribution of grid points 
in the far field which is suitable for acoustic wave propagation and the implementation of 
far field, radiation type boundary conditions. 

In order to satisfy these requirements, we make a transformation into the elliptic 
coordinates ( 77 , £) with the transformation 

$ = a* cos( 7 r 77 )cosh( 7 r£) (3.50a) 

= a* sin( 7 r? 7 )sinh( 7 r£). (3.50b) 

where a* is an arbitrary constant which will be defined later. Note that in the far field the 
elliptic coordinates reduce essentially to cylindrical coordinates, and that the $ — 'P plane 
is mapped into a semi-infinite strip in the 77 — £ plane. 

With this change of variables, the governing equation takes the following form: 


_, 2 ^ ( L)*] 
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+ A 1 J{jji 0 ip + Tx ^ + T 2 ^ + T 3 d ^ 2 + T 4 d ^ 2 + T 5 
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(3.51) 
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Here J(r],£) is the Jacobian of the transformation (3.50) which is given by 

J(r), £) = 7r 2 a* 2 [sin 2 (nr]) + sinh 2 (7r£)] , (3.52) 

and the T, coefficients are known functions of 77 and 

The airfoil boundary condition (3.48) has also been transformed through the change 
of variables (3.50). Expressing condition (3.48) in the variables (r;,£), one gets 


dip * • f \ elK °* r 
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e, 1 


(3.53) 


We now proceed to discuss the remaining boundary conditions. First, the wake bound- 
ary condition (3.49a) may be integrated so that it becomes 


A (ipe 


- i Kn ‘I’ 


) = A <p 2 + [A(if>e ! Ko *)t. e . ~ A <p 2t .Je 


ik 


*t.e. 


u o 


(3.54) 


where the subscript t.e. denotes quantities at the airfoil trailing edge. Note that, in general, 
both ip and $ are discontinuous across the wake, so that in evaluating A(f/>e - * Ao$ ), it is 
necessary to take into account the jump in both xp and <3>. This condition is imposed for 
grid points on the lower side of the wake. 

On the upper side of the wake, we impose condition (3.49c), which specifies that the 
normal velocity component of the unsteady velocity is continuous across the wake. This 
may be written 


u„ 



u ° 



(3.55) 


where and ” superscripts denote above and below the wake, respectively, and the 
derivatives in (3.55) are taken to be one-sided. 

In order to proceed further in the development of the boundary conditions (3.54) 
and (3.55) in terms of the computational coordinates (p,0) ^ necessary to first discuss 
in more detail the transformation (3.50). First, the constant a * is determined from the 
condition that the airfoil trailing edge point on the suction side, where $ = a*, T = 0, 
should map into the point r/ = 0,£ = 0, and the stagnation point, where $ = — a*,T = 0 
should map into the point 77 = 1, £ = 0. Then a* must be determined from 


a 


★ 



(3.56) 


where s denotes the arc length along the airfoil surface, t.e. denotes the airfoil trailing 
edge, and s.p. denotes the stagnation point. The steady solver FL036 is used to locate the 
stagnation point, and the integration in (3.56) is carried out using trapezoidal integration. 
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The suction surface of the airfoil, then, corresponds to the line segment on the 77 
axis between 0 and 1 , and the pressure surface corresponds to the line segment on the 
77 axis between 1 and r] t e __, where 77 < . e ._ < 2.0 (See Figure 5.) The pressure side of the 
wake then corresponds to the line segments given by {? 7 <-e ._ < 77 < 2 . 0 , £ = 0 }, and 
{77 = 2.0, { > 0}. The suction side of the wake corresponds to the positive £ axis. The 
upper boundary of the 77 — £ grid, which is given by {£ = £ max , 0 < 77 < 2 }, corresponds 
to the far field boundary in the $ — plane. 

Because of the discontinuity in $ across the wake grid line, the grid points on the 
upper and lower sides of the wake in the physical plane (the x\ — £2 plane) do not coincide. 
This presents some difficulty in the implementation of boundary conditions (3.54) and 
(3.55), inasmuch as these conditions both specify a relation that must be satisfied across 
the wake. The difficulty can be removed, however, by simply using a linear averaging of 
7/) at the two adjacent grid points to represent ip at an arbitrary point in between. Using 
such a linear averaging, then, boundary condition (3.54), which is imposed for wake grid 
points on the pressure side, becomes 





J avg 


(ipe~ iK °^) 


= Ay > 2 + [A(^e xKo *) t e - Av? 2 f .e.]e £ (3.57) 

The discontinuity in $ across the wake also leads to a shift in the location of the 
corresponding grid points in the physical plane on opposite sides of the wake. Because 
of this shift, the last several grid points on the pressure side in the physical plane extend 
past the last suction side wake grid point, so that the above averaging technique cannot 
be employed. (See Figure 5) However, because these last few points are in the far field 
where the mean flow is nearly uniform, the function ip behaves essentially as in the case 
of the thin airfoil problem, and is approximately an odd function of T. By assuming ip to 
be an odd function with respect to T in the far field, condition (3.54) becomes 
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(3.58) 


for extra far field wake grid points on the pressure side. 

For wake grid points on the suction side, the linear averaging technique can be used 
for all points, and condition (3.55) becomes 




Expressing this in 77 and £, one gets 


Uo - V2) 


avg 


(3.59) 
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(3.60) 


when the averaged values of ip lie on the rj axis. When the averaged values of ip lie on the 
right hand boundary of the rj — £ grid {r} = 2.0, £ > 0}, (3.59) becomes 
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(3.61) 


In order to complete the formulation of the boundary conditions, it is only necessary to 
specify conditions at the airfoil trailing edge and in the far field. At the trailing edge, there 
are two grid points that coincide, one corresponding to the suction side of the airfoil and one 
corresponding to the pressure side. It is therefore necessary to impose two conditions at the 
trailing edge point. At the point on the suction side, which corresponds to (??,£) = (0,0), 
the Jacobian of the coordinate transformation (3.50) vanishes. Since the Kutta condition 
requires that the velocity at the trailing edge be finite, we are led to the requirement that 



at 


( 7,0 = ( 0 , 0 ). 


(3.62) 


At the pressure side trailing edge point, we impose the condition that the unsteady 
pressure is continuous, 


p't.e.+ =P\.e.- at (??,£) = 07t.e.-,0). (3.63) 

We point out that since the Jacobian vanishes at the suction side trailing edge point, it is 
not possible to directly calculate the pressure in rj — £ coordinates at that point. However, 
it can be calculated using a Taylor series expansion from neighboring points. By using this 
approach, condition (3.63) can be satisfied. 

In presenting the far field boundary condition, we first comment that while condition 

— f 

(3.12) expresses the mathematical requirement that V (pi — * 0 at upstream infinity, this 
condition cannot be imposed throughout the far field on a boundary at a finite distance 
from the airfoil. To implement such a condition would impose a reflecting boundary 
condition which can lead to large errors in the solution. 
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To correctly model the physics of the present unsteady boundary value problem re- 
quires that the far field boundary condition be such that it allows outgoing acoustic waves 
to leave the solution domain without being reflected back into the computational grid. 
This can be accomplished, for example, by using separation of variables along with a series 
expansion for the far field solution ip and only accepting those terms in the series which 
represent outgoing waves. The difficulty with this approach, however, is that it leads to a 
matrix which requires pivoting and therefore longer solution times. In addition, since ip is 
not continuous across the wake, but the series expansion for ip is continuous everywhere, 
there is an incompatibility near the wake which can lead to a poor solution in the far field. 

An alternative to the series expansion approach is to use a Sommerfield radiation con- 
dition on the unsteady pressure. This approach avoids both the difficulties associated with 
the series expansion method and is also easier to implement. The Sommerfield radiation 
condition for the pressure is the approach used in the present work, and may be written 
in operator notation as 
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where 


$ = R cos© (3.65a) 

= i?sin 0. (3.65b) 


Neglecting terms, this reduces to 
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This condition is applied for all grid points such that 0 < 77 < 2, £ = £ max . 


(3.66) 


Numerical Method 

In the previous section we presented the transformation into computational coordi- 
nates and the development of the numerical boundary value problem. The problem to be 
solved numerically consists of the governing equation (3.51), and the boundary conditions 
(3.53), (3.57), (3.58), (3.60), (3.61), (3.62), (3.63) and (3.66). As mentioned previously, 
our basic numerical approach is to use the method of finite difference approximations, and 
then to solve the resulting linear system of equations using a matrix solver. 
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The first step in obtaining numerical solutions to equation (3.51) and its associated 
boundary conditions for a given flow configuration is to calculate the source term S and the 
coefficient functions Ai...A 5 . This requires the evaluation of Lighthill’s drift function and 
its first and second partial derivatives with respect to T, and the evaluation of the mean flow 
quantities and their partial derivatives with respect to $ and T at each interior grid point. 
The mean flow quantities are obtained through the analytical scheme presented previously 
in the present paper, and their derivatives are calculated using four point differencing. 

It should be emphasized that accurate evaluation of the source term is essential if 
accurate solutions to equation (3.51) are to be obtained, and this in turn depends largley 
upon the accurate evaluation of the drift function. If the drift function is not evaluated 
accurately, then the source term will not tend to zero in the far field and the numerical 
scheme will become unstable. One of the major advantages of using the analytical scheme 
outlined previously in this section to obtain the mean flow is that it can determine the 
mean velocities at arbitrary ($, T). This means that for fixed T, i.e. on a given streamline, 
we can determine the mean velocities for arbitrary $. Since evaluation of the drift function 
requires the integration of the expression ( ttj — jjf) with respect to $ on a fixed streamline, 
it is very easy to do the numerical integrations necessary to accurately evaluate the drift 
function A. In the actual calculations, we evaluate A at a given grid point as the sum of 
an analytically determined part and a numerically determined part. The analytical part 
comes from a far field expansion for A which is given by 

A = $ 0 - -sgn('fo) J + sgn(^o)tan~ 1 ( - $ °, ) (3.67) 

This expression can be used to accurately calculate A at some point far upstream, and 
then since A is additive, the remaining portion of the integration can be done numerically 
from the upstream location to the given grid point. The numerical integration is done 
using the trapezoid rule with variable spacing in $ to ensure accurate resolution near the 
airfoil. The first and second partial derivatives of A are approximated using four point 
differencing. 

Once S and A\...A$ have been calculated, they can be stored separately and passed 
to the subroutine which sets up the matrix equation to be solved. 

We now proceed to discuss the differencing used for the governing equation and bound- 
ary conditions. To represent equation (3.51) we use the standard nine-point, central dif- 
ference computational molecule which is second order accurate in rj and £. We assume in 
general that the spacing in each direction is nonuniform. Details of the grid spacing will be 
discussed momentarily. Each of the boundary conditions (3.53), (3.62), (3.63), and (3.66) 
are implemented using four-point, one-sided differencing which is third order accurate for 
(3.53), (3.62), and (3.63), and second order accurate for (3.66). Boundary conditions (3.60) 
and (3.61) are both implemented using three-point, one-sided differencing which is second 
order accurate. 

Obtaining a numerical solution to the finite difference equations representing the gov- 
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matrix equation whose size is equal to twice the number of grid points. There are difficul- 
ties in solving this linear system of equations because the matrix is not block tridiagonal 
and does not have a regular block structure which can be exploited. In addition, iterative 
solvers have convergence problems because the diagonal dominance of the matrix changes 
as the parameters Moo, &i, and k^ are varied. 

Because of these difficulties, a general purpose sparse matrix solver was developed 
which stores only the nonzero entries of the matrix and can solve an arbitrary sparse matrix 
equation using Gaussian elimination. The solver basically works by using an ordered list 
to represent the nonzero entries of each row in the matrix, and then inserts and deletes 
new entries in the rows of the matrix as multiples of each row are added to other rows to 
carry out the elimination process. The only requirement for the solver to work is that the 
matrix must be arranged such that it remains reasonably sparse during the elimination 
procedure. 

The sparse solver has both a pivoting and non-pivoting feature. However, as pivoting 
during the elimination process proved to be unnecessary, the pivoting feature was not 
used. By not pivoting during the elimination, it was possible to increase the storage 
efficiency of the solver and thereby solve larger systems of equations. The increased storage 
efficiency was gained by using a mapping function to map sub-blocks of the rectangular 
two-dimensional arrays containing the nonzero entries of the matrix and their associated 
column numbers into singly dimensioned arrays which contained less unused storage. By 
using this technique, the storage efficiency of the solver was increased by about 25 percent. 

The final issue to be discussed in regard to our numerical scheme is the method of grid 
determination. As reported previously , 10 ’ 11 it is not possible to use a single grid and obtain 
accurate solutions to the gust response problem for a large range of reduced frequencies. 
Rather, the unsteady grid must be determined as a function of both the mean flow Mach 
number and the reduced frequency. 

This requirement is dictated by both the accuracy of the far field boundary condition 
(3.66), and the need to adequately model the airfoil boundary condition (3.56) and the wake 
boundary conditions (3.57) and (3.58). The accuracy of the far field boundary condition 
depends on the reduced frequency k\ and free stream Mach number Moo in such a way that 
the parameter fcl 3f°° i?,, where R is the distance to the far field boundary, should remain 

'OO 

at least 0(1). This shows that the location of the outer boundary of the grid must be 
determined as a function of k\ and In addition, there should be enough grid points 

per wavelength to accurately represent the airfoil and wake boundary conditions. Due to 
the harmonic terms containing the parameters K 0 and k ly this shows that the r\ and £ 
spacing have to be determined as a function of k\ and Moo- 

The spacing in the r) direction is constant for 0 < rj < 1, and then changes slightly but 
constant again for 1 < rj < rjt. e .-, and finally constant again, but slightly different from 
the two previous intervals, for rjt.e.- < ij < 2. The spacing on 0 < rj < r] t . e .- determines 
the spacing on the airfoil surface. Normally the number of grid points in the rj direction 
varies from 40 for the low frequencies up to about 70 for reduced frequencies of 4. 


77 



An optimal spacing of the grid points in the wake (the £ direction) turns out to be 
more difficult to achieve than the spacing on the airfoil. Numerical studies of the thin 
airfoil gust response problem showed that the optimal choice of spacing was 12 uniformly- 
spaced grid points per wavelength. For the case of the thin airfoil, the wake boundary 
condition analogous to condition (3.57) is much simpler, and it is possible to choose the 
spacing of the £ grid points such that they are uniformly spaced along the waves. However, 
for the general problem of nonuniform mean flows, the waves in the wake are distorted due 
ik, f* ^ 

to the e J *t.e. u 0 term, and it is no longer possible to determine the grid point spacing 
such that they are uniformly spaced along the waves. This does not prove to be a difficulty, 
however, because in the far field the general problem of vortical flows past a lifting airfoil 
reduces essentially to a linear problem as the mean flow tends to become uniform. So we 
may determine the wake spacing as in the case of the flat plate airfoil, and in the far field 
the grid points will be nearly uniformly spaced along the waves. 

For a flat plate airfoil, the wake boundary condition which imposes the continuity of 
the pressure is given, corresponding to the transformation (3.50), by 


i4f-[a*cosh(7r^)-“*] 

Ytvake,j — Yt.e.& 00 


(3.68) 


In order to have uniformly spaced grid points along the waves in the wake, the argument 
of the exponential function should vary by equal fractional increments of 7 r. To place 12 
grid points per wavelength, we are then led to the requirement that the location of the 
j th £ grid point be determined from the relation 

Jc-t O7 r 

|K«*coshOrfc)-a*] = j~. (3.69) 

Solving for £j, we get 

+ (3 - 70) 

This is the method for determining the £ spacing in the far field. Near the airfoil the above 
procedure leads to a spacing which is too coarse to be used. So near the airfoil we use 
uniform £ spacing which remains constant at some value A£, until a point is reached such 
that the £j determined from (3.70) satisfy £j+i — £j < A£. From that point on, the spacing 
of the £ grid points is determined from (3.70). 


IV. CODE VALIDATION 

Extensive efforts were taken to validate the computer codes which were developed to 
implement the numerical solution procedure which was outlined in the previous section. 
The validation process consisted of a combination of comparing numerical results with 
known analytical solutions to the classical thin airfoil gust response problems, comparing 
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with the second order theory of Goldstein and Atassi 3 and Atassi, 4 comparing with the first 
order numerical results of Atassi and Grzedzinski, 13 and calculating solutions to limiting 
case problems i.e., as Mach number, thickness, angle of attack, or camber go to zero. 

Sample computation times for the results presented varied considerably, depending on 
the reduced frequency, Mach number, and airfoil loading. For thin, unloaded airfoils, with 
low reduced frequency gusts, typical solution times were about 20 seconds per frequency 
on the Cray X-MP at the NASA Lewis Research Center. The higher frequencies for these 
airfoils required on the order of 60 CPU seconds per frequency, with slightly higher solution 
times as the Mach number increases. For thick, symmetric, unloaded airfoils, the solution 
times range from about 40 seconds for the lower frequencies up to about 150 seconds for 
the highest frequencies. Finally, for loaded airfoils, the solution times ranged from about 
250 seconds for the low frequencies up to around 1200 seconds for the highest reduced 
frequencies. No effort was made to optimize the computational efficiency of the scheme, 
as our main purpose was to validate its accuracy. 

In the results that follow, comparisons are made for one- dimensional (transverse) 
gusts, two-dimensional (transverse and longitudinal) gusts, and fully three-dimensional 
gusts. See Figures 6, 7, and 8. 

The first step in the validation process was to compare numerical results with known 
solutions to the classical thin airfoil gust response problems. In Figures 9 and 10 we present 
comparisons between numerical and analytical results for the normalized unsteady lift for 
vortical flows past flat plate airfoils with zero thickness. The normalized unsteady lift, 
commonly called the response function, is defined by 


RL(ki,fa,Moo) 
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(4.1) 


where L' is the unsteady lift. Figure 9. a shows a comparison between numerical results and 
the Sears solution 14 for the case of a one-dimensional (transverse) gust in incompressible 
flow, and Figures 9.b and 9.c compare numerical results and analytical results obtained 
from a Possio solver for a one- dimensional gust at Mach numbers of .5 and .8. The reduced 
frequency values at which the comparisons are made range from 0 to 4.0, and are shown 
below the plot. The point on the real axis and furthest to the right corresponds to Aq = 0, 
and the other points along the curve correspond in order to the other reduced frequency 
values. Figures 10. a through 10. c compare numerical results and analytical results from a 
Possio solver for three-dimensional gusts for Mach numbers of .1, .5, and .8. The conditions 
on the gust wave number parameters are shown below the plots. As can be readily seen, 
there is excellent agreement between the numerical and analytical results. The only loss 
of accuracy occurs when both the Mach number and reduced frequency become large. 

In order to assess the accuracy of the present numerical scheme for vortical flows 
around thin airfoils in which the mean flow is not uniform, we compare with the second 
order theory of Goldstein and Atassi 3 and Atassi. 4 The results given by Atassi assume a 
zero thickness airfoil, but account for the effects of airfoil camber and angle of attack on the 
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airfoil unsteady response. In Figure 11, we compare the numerically computed response 
function with the second order theory for an incompressible flow with a two-dimensional 
gust about an airfoil with an angle of attack of two degrees and a camber ratio of .05. 
The numerically computed response function is for a 6% thick Joukowski airfoil, while 
the second order theory does not take into account the airfoil thickness. The reduced 
frequency values used for the comparison are shown below the plot. As can be seen, the 
numerical results for the 6% thick airfoil show a slightly larger lift at the low frequencies, 
but a slightly smaller lift at the higher frequencies for k\ up to about 3.0. As shown in 
the results presented in [5] and [15], this effect can be attributed entirely to the thickness 
of the airfoil, so that the agreement is very good for reduced frequencies ranging from 0 
up to about 3.0. For the frequencies higher than 3.0, it is not possible to make any firm 
conclusion on the accuracy of the numerical results. 

To validate the numerical scheme for airfoils with thickness, we compare with the 
results of Atassi and Grzedzinski. 13 In Figure 12, we show comparisons for one-, two-, 
and three-dimensional gusts for incompressible flows around a 12% thick Joukowski airfoil 
with zero degrees angle of attack and zero camber ratio. The reduced frequency values for 
the comparison range from .2 to 2.5 and are shown on the plots. We limit the comparison 
to this range of k\, since this is roughly the range of validity of the Green’s function 
approach of Atassi and Grzedzinski. The agreement between the two sets of results is 
good in general. 

The final step in the validation process was to calculate the solutions to various limiting 
case problems. The limiting case of M*, — » 0, i.e., the incompressible case, was covered 
above where the numerical results were compared to the Sears solution. We now present 
results for the limiting cases of airfoil thickness, angle of attack, and camber. 

In Figures 13 through 16, we compare numerical results for the unsteady lift and 
moment about the airfoil center of a 3% thick, symmetric Joukowski airfoil with that of 
a flat plate airfoil with zero thickness. Analogous to the response function Ri for the 
unsteady lift, we define the response function Rm for the unsteady moment by 


JZmOi^Moo) 


M' 
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where M' is the unsteady moment about the airfoil center. Figures 13 and 14 present 
results for = -1) and Figures 15 and 16 present results for = .6. For both 
Mach numbers, it is seen that the small airfoil thickness has little effect on the unsteady 
response, except for high reduced frequencies in the two-dimensional gust case, where the 
magnitude of the unsteady lift is reduced by 15 - 20 percent. It would appear from these 
results that thickness effects become more important at the higher frequencies for the case 
of the two-dimensional gust. 

Figures 17 and 18 present comparisons between results for a 12% thick symmetric 
Joukowski airfoil at zero degrees angle of attack and one degree angle of attack. All plots 
are for a free stream Mach number of .1. As in the case of airfoil thickness, the strongest 
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effect is seen in the two-dimensional gust case. However, here there is a significant effect 
both for the low and high reduced frequencies. At the low frequency end, the effect is 
primarily a reduction in the magnitude of the unsteady lift, while at the high end it is 
primarily a change in phase of the unsteady lift. We also point out that, in agreement 
with the theoretical results of Atassi, 4 for the transverse gust case in which the gust has 
only an upwash component, the steady loading on the airfoil has virtually no effect on the 
unsteady lift. 

Finally, in Figures 19 through 20, we compare results for a 12% thick Joukowski airfoil 
with no camber with a 12% thick Joukowski airfoil with a camber ratio of .02. The free 
stream Mach number for all plots is .1, and the angle of attack is zero degrees. The effect 
here is exactly analogous to the effect of small angle of attack, except that it is stronger 
in this case due to the increased steady loading on the airfoil. For the one degree angle of 
attack airfoil the steady lift coefficient was .12, while for the 2% cambered airfoil the steady 
lift coefficient was .27. In each case, the reduction of the quasi-steady lift (Aq = 0 = Aq) for 
the two-dimensional gust is directly proportional to the steady loading on the airfoil with 
a proportionality constant of 0.26. Using the theoretical results of Atassi reported in [4], 
it can be shown that for zero thickness airfoils in a two-dimensional gust in incompressible 
flow, the reduction in the quasi-steady lift for airfoils with small camber and angle of 
attack is proportional to the steady lift coefficient with a proportionality constant of ^ 

= .23. The difference between the numerical and theoretical values of the proportionality 
constant can be accounted for by the fact that the theoretical result does not account for 
the thickness of the airfoil, and also assumes a parabolic profile for the airfoil camber line. 

In Reference 15 a more general study of how the two-dimensional quasi-steady lift 
changes in response to mean airfoil loading showed that for airfoils with heavier loading 
and Mach number effects the value of the proportionality constant was about .25. We thus 
conclude that the reduction in the two-dimensional quasi-steady lift due to mean airfoil 
loading is roughly proportional to the steady lift coefficient of the loaded airfoil with a 
proportionality constant of .25. 

Before concluding the discussion of numerical results, the authors would like to em- 
phasize the significance of the method of grid determination which was outlined in the 
previous section. In Figure 21 we present numerical results which demonstrate the kinds 
of errors that can occur as a result of using an inappropriate grid. The results shown in 
these figures were generated without determining the grid as a function of the reduced 
frequency. For each case shown, the same grid was used for all frequencies in the calcula- 
tion. The grid used for each Mach number was the one normally used only for the highest 
reduced frequency. By using the grid for the highest frequency, it was assured that there 
would be sufficient grid resolution to resolve the waves for the lower frequencies. But as can 
be seen, the agreement is not nearly as good as when the grid is determined as a function 
of both the Mach number and reduced frequency. These results show the effect of keeping 
the outer boundary fixed, and not varying it with the reduced frequency in order to ensure 
that the representation of the far field boundary condition is sufficiently accurate. 
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The results in Figure 22 show the kinds of errors that can occur when the grid points 
are not suitably spaced in the far field. The grids used for the results in Figure 22 used 
an exponentially decreasing spacing which was varied to ensure that there were enough 
grid points per wavelength to adequately model the wavelike structure of the solution. In 
addition, the location of the far field boundary was also varied to ensure that the far field 
boundary condition would be sufficiently accurate. However, as the results show, there 
are large errors in the response function curves. This is due to the fact that exponential 
spacing is not suitable for this kind of wave propagation problem. 


V. CONCLUSION 

In the present paper the authors have presented a finite-difference, frequency- domain 
numerical scheme for the solution of unsteady, subsonic vortical flows around lifting airfoils. 
The present method is an alternative to the potential and primitive variable methods, and 
due to its inherent efficiency and ability to correctly handle the convection and distortion of 
the upstream vorticity by the nonuniform mean flow, represents an ideal solution method 
for unsteady aerodynamic flow fields. 

The computer codes that have been developed to implement the present numerical 
approach have been validated through extensive comparisons with known solutions to un- 
steady vortical flow problems and through the calculation of various limiting case problems. 
The numerical results have shown that our numerical scheme can calculate with very good 
accuracy the solutions to a large variety of unsteady vortical flow problems. We conclude 
that for symmetric airfoils and loaded airfoils with small mean loading, the present scheme 
is very accurate for reduced frequencies ranging from 0 up to at least 4.0, for both incom- 
pressible and compressible flows. For more heavily loaded airfoils, the scheme has similar 
accuracy for reduced frequencies ranging from 0 up to around 3.0. For reduced frequencies 
above 3.0, it is not possible to make any firm conclusion on the accuracy of the present 
numerical scheme for airfoils with substantial mean loading. This is an area that requires 
further study. 

Among the most important features of our numerical approach are the transformation 
of the independent variables into elliptic coordinates, the method of determining the un- 
steady grid as a function of the Mach number and reduced frequency, the far field radiation 
condition for the unsteady pressure, the general purpose direct sparse matrix solver, and 
the formulation and method of evaluation of the source term. 

Finally, the authors are in the process of extending the present linearized unsteady 
aerodynamic analysis and numerical solution procedure to include transonic flows. Details 
will be presented in a future paper. 
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MACH NUMBER 



Figure 1. Comparison of Mach number at the airfoil surface between 

FL036 and the present analytical scheme. = .6, a = 0°, 

camber = 0, thickness = .12. 
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Figure 2. Comparison of Mach number at the airfoil surface between 
FL036 and the present analytical scheme. M ^ — .6, a = z 
camber = 0, thickness = .12. 
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Figure 3. Comparison of Mach number along a FL036 grid line which 
wraps around the airfoil. — .6, a = 0°, camber = 0, 

thickness = .12. 
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Figure 4. Comparison of Mach number along a FL036 grid line which 
wraps around the airfoil. Moo — .6, a — 2°, camber = 0, 
thickness = .12. 
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Figure 5. Computational and physical grids. 
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Figure 6. Airfoil in a transverse gust. 




Figure 7. Airfoil in a transverse and longitudinal gust 




Figure 8 . Airfoil in a three-dimensional gust. 
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REAL LIFT 


Figure 9. Comparison between the numerically computed unsteady lift 
and analytical results for a flat plate airfoil in a transverse 
gust at (a) M = 0.1, (b) M = 0.5, and (c) M = 0.8 
kx = 0.0, 0.007, 0.027, 0.062, 0.110, 0.172, 0.24S, 0.338, 0.442, 0.561, 
0.694, 0.S42, 1.01, 1.18, 1.38, 1.59, 1.82, 2.07, 2.33, 2.62, 2.93, 
3.26, 3.62, 4.01 
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Figure 10. Comparison between the numerically computed unsteady lift 
and analytical results for a flat plate airfoil in a three- 
dimensional gust at (a) M = 0.1, (b) M — 0.5, and (c) M = 0.8 
k 3 = 0.442, |a| = 1, = — j, k\ = k 2 , a ■ k — 0, > 0 

fci = 0.0, 0.007, 0.027, 0.062, 0.110, 0.172, 0.248, 0.338, 0.442, 0.561, 
0.694, 0.842, 1.01, 1.18, 1.38, 1.59, 1.82, 2.07, 2.33, 2.62, 2.93, 
3.26, 3.62, 4.01 
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Imaginary Lift 



- 0.2 0.0 0.2 0.4 0.6 


Real Lift 

Figure 11. Comparison between the numerically computed unsteady lift 
and the second order theory for an airfoil in a transverse 
and longitudinal gust. The second order theory does not account for 
the thickness of the airfoil. Moo = .1, <* = 2°, camber = .05, 
thickness ratio = .06. —a\ = a 2 — .7071, ki = ^ 2,03 — kz = 0 
k\ = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Real Lift 

Figure 12. Comparison between the numerically computed unsteady lift and the 
first order results of Atassi and Grzedzinski for a (a) transverse gust, 
(b) transverse and longitudinal gust with —a\ = a 2 = .7071, k\ = & 2 , 
a 3 = k 3 = 0, and (c) a three-dimensional gust with k 3 = 0.4, \a\ = 1, 

= _I ? kx = k 2 , a ■ k = 0, a 2 > 0. = .1, a = 0, camber = 0, 

thickness ratio = .12. Results shown are the complex conjugate values 
of the unsteady lift. 

k x = 0.2, 0.3, 0.45, 0.6, 0.8, 1.0, 1.3, 1.6, 2.0, 2.5 
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o Thickness = .03 
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Figure 13. Comparison between the unsteady lift of a flat plate airfoil and a 
3 percent thick Joukowski airfoil in a (a) transverse gust, 

(b) transverse and longitudinal gust with —a% — a 2 — .7071, k\ = } 
a 3 = = 0, and (c) a three-dimensional gust with = 0.4, |a| = ' 

^ = -7, hi - k 2 , a - k = 0, a 2 > 0. = .1, a — 0, camber = 0 

= 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure 14. 


Comparison between the unsteady moment of a flat plate airfoil and a 
3 percent thick Joukowski airfoil in a (a) transverse gust, 

(b) transverse and longitudinal gust with —ai = a 2 — .7071, k\ = h?, 
a 3 — kz =0, and (c) a three-dimensional gust with £3 = 0.4, |a| = 1, 
^ fci = k 2 , a ■ k — 0 , a 2 > 0 . Mqo = - 1 , a = 0 , camber = 0 . 

fcj = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure 15. Comparison between the unsteady lift of a flat plate airfoil and a 
3 percent thick Joukowski airfoil in a (a) transverse gust, 

(b) transverse and longitudinal gust with —a i = 02 = .7071, k\ - 
a 3 = k 3 = 0, and (c) a three-dimensional gust with k 3 = 0.4, \a\ - 
si. - -1 ki = ko, a-ic = 0, a o > 0. = .6, a = 0, camber = 

h = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure 16. Comparison between the unsteady moment of a flat plate airfoil and a 
3 percent thick Joukowski airfoil in a (a) transverse gust, 

(b) transverse and longitudinal gust with — a\ = a 2 = .7071, k\ = fc 2 , 
a 3 = fc 3 = 0, and (c) a three-dimensional gust with k 3 = 0.4, |a| = 1, 
= - 7, ki - k 2 , a ■ k = 0, a 2 > 0. M <*> = .6, a — 0, camber = 0. 
ki = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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REAL LIFT 

Figure IT. Comparison between the unsteady lift of a 12 percent thick, Joukowski 
airfoil at 0° angle of attack and 1° angle of attack for a (a) transverse 
gust, (b) transverse and longitudinal gust with —a x = a 2 = .7071, k x - 
a 3 = k 3 = 0, and (c) a three-dimensional gust with k 3 = 0.4, |a| = 1, 
^ = -7, k x = k 2 , a-k — 0, a 2 > 0. = .1, camber = 0, 

thickness ratio = .12. 

ki = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.S, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure IS. Comparison between the unsteady moment of a 12 percent thick, Joukowski 
airfoil at 0 ° angle of attack and 1° angle of attack for a (a) transverse 
gust, (b) transverse and longitudinal gust with —a i = <22 = .7071, k\ — ko, 
a 3 = I’ 3 = 0, and (c) a three-dimensional gust with £3 = 0.4, |a[ = 1, 
tti = — 1 , ki = k 2 , a ■ k = 0 , a 2 > 0 . M ^ — . 1 , camber = 0 , 
thickness ratio = . 12 . 

^ = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure 19. Comparison between the unsteady lift of an uncambered Joukowski 
airfoil and an airfoil with camber ratio of .02 for a (a) transverse 
gust, (b) transverse and longitudinal gust with — a\ — a 2 = .7071, ki = k 2 , 
0-3 = k 3 = 0, and (c) a three-dimensional gust with kz = 0.4, ja| = 1, 

^ — -j, ki = k 2 , a • k — 0, a.2 > 0. Moo — -1) a — 0°i 
thickness ratio = .12. 


h = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 
0.8, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure 20. Comparison between the unsteady moment of an uncambered Joukowski 
airfoil and an airfoil with camber ratio of .02 for a (a) transverse 
gust, (b) transverse and longitudinal gust with — a\ = <12 — .7071, k\ — ko, 
a 3 — = 0, and (c) a three-dimensional gust with k 3 = 0.4, |a| = 1, 

^ = — j, ki = k 2 , a ■ k = 0, a 2 > 0. M ^ = .1, a; = 0°, 
thickness ratio = .12. 

ki = 0.0, 0.01, 0.03, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 

0.S, 1.0, 1.3, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0 
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Figure 21. Comparison between numerical results generated on a fixed grid 

and analytical results for a flat plate airfoil in a transverse gust at 
(a) M — 0.1, (b) M — 0.5, and (c) M = 0.8. 
ky = 0.0, 0.007, 0.027, 0.062, 0.110, 0.172, 0.248, 0.338, 0.442, 
0.561, 0.694, 0.842, 1.01, 1.18, 1.38, 1.59, 1.82, 2.07, 2.33, 
2.62, 2.93, 3.26, 3.62, 4.01 


103 



0 s QQ 


• Moo = .2 

■ Moo = .5 


g„& 4 r 



REAL LIFT 


Figure 22. Numerical results for the unsteady lift of a 12 percent thick, 

symmetric Joukowski airfoil in a transverse gust. Results show grid 
dependent errors due to exponentially decreasing spacing in the 
£ direction which was used in place of condition (3.70). 

k x = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 

1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 

2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0 
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